#ifndef __fe_evaluators_hh__
#define __fe_evaluators_hh__

#include "femlib_config.hh"
#include "point.hh"
#include "fem_types.hh"
#include "fem_traits.hh"
#include "femlib_exceptions.hh"
//#include "fe.hh"

#if defined (FADBADPP_GRADIENTS) || defined (FADBADPP_HESSIANS)
#include "fadiff.hh"
#include "badiff.hh"
#include "tadiff.hh"
#endif


template<class T, class Imp>
struct BasisEvaluator
{
  public:
    typedef T BasisDataType;
    typedef Imp FiniteElementType;
    
    /** virtual function which evaluates the values of all
      shape functions at a given point ip */
    BasisDataType operator() ( unsigned int ib, const Point<SPATIAL_DIMENSION, BasisDataType >& pt)
    {
      return FiniteElementType::template evaluate_basis< BasisDataType > (ib, pt);
    }

    /** virtual function which evaluates the values of all
         shape functions at all point ips */
    void operator() ( const PointList<SPATIAL_DIMENSION, BasisDataType >& ips, std::vector< BasisDataType >& basis_vals )
    {
      FiniteElementType::template evaluate_basis_all< BasisDataType > (ips, basis_vals);
    }
    
};


#ifdef FADBADPP_GRADIENTS

  template<typename T>
  struct FEGradientTraits
  {
    typedef fadbad::F<T> FADGradient;
    enum { FADGradient_Type = fastfem::FORWARD } ;
  };

  template<class Imp, class T, class U=typename FEGradientTraits<T>::FADGradient>
  struct GradientEvaluator
  {
    public:
      typedef T BasisDataType;
      typedef U GradientDataType;
      typedef Imp FiniteElementType;
      
      /** virtual function which evaluates the values of all
        shape functions at a given point ip */
      GradientDataType evaluate_basis ( unsigned int ib, const Point<SPATIAL_DIMENSION, GradientDataType >& pt)
      {
        return FiniteElementType::template evaluate_basis< GradientDataType > (ib, pt);
      }

      /** virtual function which evaluates the values of all
           shape functions at all point ips */
      void evaluate_basis_all ( const PointList<SPATIAL_DIMENSION, GradientDataType >& ips, std::vector< GradientDataType >& basis_vals )
      {
        FiniteElementType::template evaluate_basis_all< GradientDataType > (ips, basis_vals);
      }
      
      void operator() ( unsigned int ib, const Point<SPATIAL_DIMENSION, BasisDataType >& pt, Point<SPATIAL_DIMENSION, BasisDataType >& grad)
      {
        Point<SPATIAL_DIMENSION, GradientDataType > pp(pt);

        // Differentiate with respect to each primary variable
        for (unsigned int id=0; id < SPATIAL_DIMENSION; ++id)
          pp[id].diff(id, SPATIAL_DIMENSION);
        
        // Evaluate function and derivatives
        GradientDataType f = this->evaluate_basis(ib, pp);
        //GradientDataType f = FiniteElementType::template evaluate_basis< GradientDataType >(ib, pp);
        // f.diff(0, 1) ;
        
        // Value of basis function at the requested point
        //double fval=f.x();
        
        // store the value of derivatives: df/dx into input variable 'grad'
        for (unsigned int id=0; id < SPATIAL_DIMENSION; ++id)
          grad[id] = f.d(id);

        return ;
      }
      
      void operator() ( unsigned int ib, const Point<SPATIAL_DIMENSION, BasisDataType >& pt, BasisDataType& basis_val, Point<SPATIAL_DIMENSION, BasisDataType >& grad)
      {
        Point<SPATIAL_DIMENSION, GradientDataType > pp(pt);

        // Differentiate with respect to each primary variable
        for (unsigned int id=0; id < SPATIAL_DIMENSION; ++id)
          pp[id].diff(id, SPATIAL_DIMENSION);
        
        // Evaluate function and derivatives
        GradientDataType f = this->evaluate_basis(ib, pp);
        //GradientDataType f = FiniteElementType::template evaluate_basis< GradientDataType >(ib, pp);
        // f.diff(0, 1) ;
        
        // Value of basis function at the requested point
        basis_val = f.x();
        
        // store the value of derivatives: df/dx into input variable 'grad'
        for (unsigned int id=0; id < SPATIAL_DIMENSION; ++id)
          grad[id] = f.d(id);

        return ;
      }
      
      void operator() ( const PointList<SPATIAL_DIMENSION, BasisDataType>& ips, std::vector<Point<SPATIAL_DIMENSION, BasisDataType> >& grad_vals )
      {
        grad_vals.resize( ips.n_points * FEMTraits<Imp>::N_Nodes ) ;

        for (unsigned int in = 0, offset=0; in < ips.n_points; ++in, offset+=FEMTraits<Imp>::N_Nodes )
        {
          Point<SPATIAL_DIMENSION, GradientDataType > pp(ips.get_point(in));

          // Differentiate with respect to each primary variable
          for (unsigned int id=0; id < SPATIAL_DIMENSION; ++id)
            pp[id].diff(id, SPATIAL_DIMENSION);

          for (unsigned int ib = 0; ib < FEMTraits<Imp>::N_Nodes; ++ib)
          {
            // Evaluate function and derivatives
            GradientDataType f = this->evaluate_basis(ib, pp);
            //GradientDataType f = Imp::evaluate_basis<GradientDataType>(ib, pp);
            // f.diff(0, 1) ;
            
            // Value of basis function at the requested point
            //double fval=f.x();
            
            // store the value of derivatives: df/dx into input variable 'grad'
            for (unsigned int id=0; id < SPATIAL_DIMENSION; ++id)
              grad_vals[offset+ib][id] = f.d(id);
          }
        }

        return ;
      }
      
      void operator() ( const PointList<SPATIAL_DIMENSION, BasisDataType>& ips, 
                        std::vector<BasisDataType>& basis_vals,
                        std::vector<Point<SPATIAL_DIMENSION, BasisDataType> >& grad_vals)
      {
        basis_vals.resize( ips.n_points * FEMTraits<Imp>::N_Nodes ) ;
        grad_vals.resize( ips.n_points * FEMTraits<Imp>::N_Nodes ) ;

        for (unsigned int in = 0, offset=0; in < ips.n_points; ++in, offset+=FEMTraits<Imp>::N_Nodes )
        {
          Point<SPATIAL_DIMENSION, GradientDataType > pp(ips.get_point(in));

          // Differentiate with respect to each primary variable
          for (unsigned int id=0; id < SPATIAL_DIMENSION; ++id)
            pp[id].diff(id, SPATIAL_DIMENSION);

          for (unsigned int ib = 0; ib < FEMTraits<Imp>::N_Nodes; ++ib)
          {
            // Evaluate function and derivatives
            GradientDataType f = this->evaluate_basis(ib, pp);
            //GradientDataType f = Imp::evaluate_basis<GradientDataType>(ib, pp);
            // f.diff(0, 1) ;
            
            // Value of basis function at the requested point
            basis_vals[offset+ib] = f.x();
            
            // store the value of derivatives: df/dx into input variable 'grad'
            for (unsigned int id=0; id < SPATIAL_DIMENSION; ++id)
              grad_vals[offset+ib][id] = f.d(id);
          }
        }

        return ;
      }
  };
  
#endif


#ifdef FADBADPP_HESSIANS

  template<typename T>
  struct FEHessianTraits
  {
    typedef fadbad::B<fadbad::F<T> > FADHessian;
    enum { FADHessian_Type = fastfem::BACKWARD_FORWARD } ;
  };

  template<class Imp, class T, class U=typename FEHessianTraits<T>::FADHessian>
  struct HessianEvaluator
  {
    public:
      typedef T BasisDataType;
      typedef U HessianDataType;
      typedef Imp FiniteElementType;
      
      /** virtual function which evaluates the values of all
        shape functions at a given point ip */
      HessianDataType evaluate_basis ( unsigned int ib, 
                                        const Point<SPATIAL_DIMENSION, 
                                        HessianDataType>& pt)
      {
        return FiniteElementType::evaluate_basis(ib, pt);
      }

      /** virtual function which evaluates the values of all
           shape functions at all point ips */
      void evaluate_basis_all ( const PointList<SPATIAL_DIMENSION, 
                                  HessianDataType>& ips, 
                                  std::vector<HessianDataType>& basis_vals )
      {
        FiniteElementType::evaluate_basis_all(ips, basis_vals);
      }
      
      void operator() ( unsigned int ib, 
                         const Point<SPATIAL_DIMENSION, BasisDataType>& px, 
                         PointList<SPATIAL_DIMENSION, BasisDataType>& hess)
      {
        femlib_assert(ib < FEMTraits<Imp>::N_Nodes, "Invalid basis index request");
        
        Point<SPATIAL_DIMENSION, HessianDataType> pp(px);

        // Differentiate with respect to each primary variable
        for (unsigned int id=0; id < SPATIAL_DIMENSION; ++id)
          pp[id].x().diff(id, SPATIAL_DIMENSION);   // (index id of SPATIAL_DIMENSION)
                
        // Evaluate function and derivatives
        HessianDataType f = evaluate_basis(ib, pp);
        f.diff(0, 1) ;
        
        // Value of basis function at the requested point
        // double fval=f.x().x();
        // double dfdx=x.d(0).x();    // Value of df/dx
        // double dfdy=y.d(0).x();    // Value of df/dy
        // double dfdxdx=x.d(0).d(0); // Value of df/dx
        // double dfdxdy=x.d(0).d(1); // Value of df/dxdy
        // double dfdydx=y.d(0).d(0); // Value of df/dydx
        // double dfdydy=y.d(0).d(1); // Value of df/dy
        
        // store the value of derivatives: df/dx into input variable 'grad'
        for (unsigned int id=0; id < SPATIAL_DIMENSION; ++id)
        {
          T* pval = hess.get(id);
          
          for (unsigned int jd=0; jd < SPATIAL_DIMENSION; ++jd)
            pval[jd] = pp[id].d(id).d(jd);
        }
        
        return ;
      }
      
      void operator() ( unsigned int ib, 
                         const Point<SPATIAL_DIMENSION, BasisDataType>& px,
                         BasisDataType& basis_val, 
                         Point<SPATIAL_DIMENSION, BasisDataType>& grad,
                         PointList<SPATIAL_DIMENSION, BasisDataType>& hess)
      {
        femlib_assert(ib < FEMTraits<Imp>::N_Nodes, "Invalid basis index request");
        
        Point<SPATIAL_DIMENSION, HessianDataType> pp(px);

        // Differentiate with respect to each primary variable
        for (unsigned int id=0; id < SPATIAL_DIMENSION; ++id)
          pp[id].x().diff(id, SPATIAL_DIMENSION);   // (index id of SPATIAL_DIMENSION)
                
        // Evaluate function and derivatives
        HessianDataType f = evaluate_basis(ib, pp);
        f.diff(0, 1) ;
        
        // Value of basis function at the requested point
        // double fval=f.x().x();
        // double dfdx=x.d(0).x();    // Value of df/dx
        // double dfdy=y.d(0).x();    // Value of df/dy
        // double dfdxdx=x.d(0).d(0); // Value of df/dx
        // double dfdxdy=x.d(0).d(1); // Value of df/dxdy
        // double dfdydx=y.d(0).d(0); // Value of df/dydx
        // double dfdydy=y.d(0).d(1); // Value of df/dy
        
        basis_val = f.x().x();
        for (unsigned int id=0; id < SPATIAL_DIMENSION; ++id)
          grad[0] = pp[id].d(0).x();   // (derivative index id of SPATIAL_DIMENSION)
        
        // store the value of derivatives: df/dx into input variable 'grad'
        for (unsigned int id=0; id < SPATIAL_DIMENSION; ++id)
        {
          T* pval = hess.get(id);
          
          for (unsigned int jd=0; jd < SPATIAL_DIMENSION; ++jd)
            pval[jd] = pp[id].d(id).d(jd);
        }
        
        return ;
      }
      
      void operator() ( const PointList<SPATIAL_DIMENSION, BasisDataType>& ips, 
                        std::vector<PointList<SPATIAL_DIMENSION, BasisDataType> >& hess_vals )
      {
        hess_vals.resize( ips.n_points * FEMTraits<Imp>::N_Nodes ) ;
        for (unsigned int in = 0, offset=0; in < ips.n_points; ++in, offset+=FEMTraits<Imp>::N_Nodes )
        {
          Point<SPATIAL_DIMENSION, HessianDataType > pp(ips.get_point(in));

          // Differentiate with respect to each primary variable
          for (unsigned int id=0; id < SPATIAL_DIMENSION; ++id)
            pp[id].x().diff(id, SPATIAL_DIMENSION);
       
          
          for (unsigned int ib = 0; ib < FEMTraits<Imp>::N_Nodes; ++ib)
          {
            // Evaluate function and derivatives
            HessianDataType f = evaluate_basis(ib, pp);
            f.diff(0, 1) ;

            PointList<SPATIAL_DIMENSION, BasisDataType>& hess = hess_vals[offset+ib] ;

            // Value of basis function at the requested point
            //double fval=f.x();
            
            // store the value of derivatives: df/dx into input variable 'grad'
            for (unsigned int id=0; id < SPATIAL_DIMENSION; ++id)
            {
              T* pval = hess.get(id);
              
              for (unsigned int jd=0; jd < SPATIAL_DIMENSION; ++jd)
                pval[jd] = pp[id].d(id).d(jd);
            }
          }
        }

        return ;
      }
      
      void operator() ( const PointList<SPATIAL_DIMENSION, BasisDataType>& ips, 
                        std::vector<BasisDataType>& basis_vals,
                        std::vector<Point<SPATIAL_DIMENSION, BasisDataType> >& grad_vals,
                        std::vector<PointList<SPATIAL_DIMENSION, BasisDataType> >& hess_vals )
      {
        basis_vals.resize( ips.n_points * FEMTraits<Imp>::N_Nodes ) ;
        grad_vals.resize( ips.n_points * FEMTraits<Imp>::N_Nodes ) ;
        hess_vals.resize( ips.n_points * FEMTraits<Imp>::N_Nodes ) ;
        for (unsigned int in = 0, offset=0; in < ips.n_points; ++in, offset+=FEMTraits<Imp>::N_Nodes )
        {
          Point<SPATIAL_DIMENSION, HessianDataType > pp(ips.get_point(in));

          // Differentiate with respect to each primary variable
          for (unsigned int id=0; id < SPATIAL_DIMENSION; ++id)
            pp[id].x().diff(id, SPATIAL_DIMENSION);
       
          for (unsigned int ib = 0; ib < FEMTraits<Imp>::N_Nodes; ++ib)
          {
            // Evaluate function and derivatives
            HessianDataType f = evaluate_basis(ib, pp);
            f.diff(0, 1) ;

            PointList<SPATIAL_DIMENSION, BasisDataType>& hess = hess_vals[offset+ib] ;

            // Value of basis function at the requested point
            //double fval=f.x();
            
            basis_vals[offset+ib] = f.x().x();
            for (unsigned int id=0; id < SPATIAL_DIMENSION; ++id)
              grad_vals[offset+ib][id] = pp[id].d(0).x();   // (derivative index id of SPATIAL_DIMENSION)
            
            // store the value of derivatives: df/dx into input variable 'grad'
            for (unsigned int id=0; id < SPATIAL_DIMENSION; ++id)
            {
              T* pval = hess.get(id);
              
              for (unsigned int jd=0; jd < SPATIAL_DIMENSION; ++jd)
                pval[jd] = pp[id].d(id).d(jd);
            }
          }
        }

        return ;
      }
      
  };
  
#endif



#endif // __fe_evaluators_hh__

